
#include "stdio.h"
#include "assert.h"
#include "string.h"
#include "mpfr.h"
#include "mpfi.h"

#include <caml/mlvalues.h>
#include <caml/memory.h>
#include <caml/alloc.h>
#include <caml/custom.h>

// see pool.graffle for memory pool documentation
// see ideal_pool.ml for an idealized ML implementation  
// see pool.ml for what we intend to implement in C for
// reasons of efficiency, interaction with the garbage collector, etc.

/* -------------------------------------------------------------------------- */
/*  Debug                                                                     */
/* -------------------------------------------------------------------------- */

int DEBUG = 1;

void debug_print(char* s){
  if(DEBUG){
    printf("%s\n",s);
    fflush(stdout);
  }
  return;
}

void error(char* s){
  if(DEBUG){
    printf("%s\n",s);
    fflush(stdout);
  }
  return;
}

value debug_(value v){
  CAMLparam1(v);
  DEBUG = Bool_val(v);
  CAMLreturn(Val_unit);
}

/* -------------------------------------------------------------------------- */
/*  Memory pool                                                               */
/* -------------------------------------------------------------------------- */

// A pool_elem is an array that remembers how many live/active elements it
// has.  When this number reaches 0, the array is freed.
// The memory pool is a collection of 'pool_elem's
typedef struct {
  mpfi_t* arr;
  int active; // number of active pool elements
} pool_elem;

// A location is the coordinates of a particular interval
typedef struct { 
  pool_elem* pool_elem;
  int index;
} location;


// The size of a pool_elem array
int DEFAULT_SIZE = 0;

// A global pointer to the next unused interval
// An invariant of this program is that after 
// 'mpfi_initialize' is called, NEXT always "points to" a 
// fresh, allocated interval.
location* NEXT = NULL;

// the total number of active locations
// should be a multiple of DEFAULT_SIZE
int ALLOCATED_LOCATIONS = 0;

// the total number of used locations
int USED_LOCATIONS = 0;

// See section 18.9.2 of the OCaml manual for a description of
// USED and MAX
int GC_USED = 0;  
int GC_MAX = 1;   

// create a new pool_elem when the current store is used up
pool_elem* new_pool_elem(void){
  pool_elem* p = (pool_elem*) malloc(sizeof(pool_elem));
  p->arr = (mpfi_t*) malloc(DEFAULT_SIZE * sizeof(mpfi_t));
  int i;
  for(i=0;i<DEFAULT_SIZE;i++){
    mpfi_init(p->arr[i]);
  }
  p->active = DEFAULT_SIZE;
  ALLOCATED_LOCATIONS += DEFAULT_SIZE;
  return p;
}

// initialize NEXT and DEFAULT_SIZE
void mpfi_initialize(int default_size){
  DEFAULT_SIZE = default_size;
  pool_elem* p = new_pool_elem();
  NEXT = (location*) malloc(sizeof(location));
  NEXT->pool_elem = p;
  NEXT->index = 0;
  return;
}

// Fill in the fields of loc with a valid new interval,
// and "increment" NEXT
void next_location(location* loc){
  loc->pool_elem = NEXT->pool_elem;
  loc->index = NEXT->index;
  if(NEXT->index < DEFAULT_SIZE - 1){
    NEXT->index++;
  } 
  else{
    NEXT->pool_elem = new_pool_elem();
    NEXT->index = 0;
  }
  USED_LOCATIONS++;
  return;
}

// the define macro gives %10 speed increase
#define value_to_location(v) (Data_custom_val(v))

/*
location* value_to_location(value v){
  location* loc = (location*) Data_custom_val(v);
  return loc;
}
*/

// Free the space allocted to an interval
// If the corresponding pool_elem becomes empty,
// free the array
void free_location(location* loc){
  loc->pool_elem->active--;
  if(loc->pool_elem->active == 0){
    int i;
    for(i=0;i<DEFAULT_SIZE;i++){
      mpfi_clear(loc->pool_elem->arr[i]);
    }
    free(loc->pool_elem->arr);
    free(loc->pool_elem);
    ALLOCATED_LOCATIONS -= DEFAULT_SIZE;
  }
  USED_LOCATIONS --;
  return;
}

void finalize_value(value v){
  location *loc = value_to_location(v);
  free_location(loc);
  return;
} 

static struct custom_operations mpfi_ops = {
  "mpfi",
  &finalize_value,
  custom_compare_default,
  custom_hash_default,
  custom_serialize_default,
  custom_deserialize_default
};

// create a new caml value that stores a location
value new_value(){
  value v = caml_alloc_custom(&mpfi_ops,sizeof(location),GC_USED,GC_MAX); 
  location* loc = value_to_location(v);
  next_location(loc);
  return v;
}

// macro is faster
#define location_to_mpfi(loc) (&(loc->pool_elem->arr[loc->index]));

/*
mpfi_t* location_to_mpfi(location* loc){
  mpfi_t* m = &(loc->pool_elem->arr[loc->index]);
  return m;
}
*/

/* -------------------------------------------------------------------------- */
/*  Functions                                                                 */
/* -------------------------------------------------------------------------- */
  
value used_locations_(value unit){
  CAMLparam1(unit);
  CAMLreturn(Val_int(USED_LOCATIONS));
}

value allocated_locations_(value unit){
  CAMLparam1(unit);
  CAMLreturn(Val_int(ALLOCATED_LOCATIONS));
}

value mpfi_initialize_(value size){
  CAMLparam1(size);
  mpfi_initialize(Int_val(size));
  CAMLreturn(Val_unit);
}

value set_gc_used_(value n){
  CAMLparam1(n);
  GC_USED = Int_val(n);
  CAMLreturn(Val_unit);
}

value set_gc_max_(value n){
  CAMLparam1(n);
  GC_MAX = Int_val(n);
  CAMLreturn(Val_unit);
}
  
#define BASE10 10
#define MAX_IVAL_LENGTH 200
#define MAX_EXP_LENGTH 12 /*  overestimate.  max exp is 2^11 for 64 bit floats */

value of_string_(value s) {
  CAMLparam1(s);
  CAMLlocal1(x);
  x = new_value();
  location* loc_x;
  loc_x = value_to_location(x);
  mpfi_t* m_x;
  m_x = location_to_mpfi(loc_x);  
  char* cs = String_val(s);
  mpfi_set_str(*m_x,cs,BASE10);
  CAMLreturn(x);
}

value to_string_(value x) {
  CAMLparam1(x);
  CAMLlocal1(s);
  location *loc_x;
  loc_x = value_to_location(x);
  mpfi_t* m_x;
  m_x = location_to_mpfi(loc_x);
  /* left endpoint */
  mp_exp_t exp1,exp2;
  mpfr_t x1,x2;
  mpfr_init(x1);
  mpfr_init(x2);
  mpfi_get_left(x1,*m_x);
  mpfi_get_right(x2,*m_x);
  char* cs1 = mpfr_get_str(NULL,&exp1,BASE10,0,x1,GMP_RNDN);
  char* cs2 = mpfr_get_str(NULL,&exp2,BASE10,0,x2,GMP_RNDU);
  // printf("left: %s\nright: %s\n",s1,s2);
  char exps1[MAX_EXP_LENGTH] = "";
  char exps2[MAX_EXP_LENGTH] = "";
  if(!(mpfr_nan_p(x1) || mpfr_inf_p(x1)))
    sprintf(exps1,"E%d",(int)exp1); // cast is safe, as mp_exp_t = int
  if(!(mpfr_nan_p(x2) || mpfr_inf_p(x2)))
    sprintf(exps2,"E%d",(int)exp2); // cast is safe, as mp_exp_t = int
  //printf("left: %s\nright: %s\n",exps1,exps2);
  char rets[MAX_IVAL_LENGTH] = "[";
  strcat(rets,cs1);
  strcat(rets,exps1);
  strcat(rets,",");
  strcat(rets,cs2);
  strcat(rets,exps2);
  strcat(rets,"]");
  s = caml_copy_string(rets);
  mpfr_free_str(cs1);
  mpfr_free_str(cs2);
  mpfr_clear(x1);
  mpfr_clear(x2);
  CAMLreturn(s);
}

/* -------------------------------------------------------------------------- */
/*  Binary Operators                                                          */
/* -------------------------------------------------------------------------- */

// mpfi -> mpfi -> mpfi

#define BINOP void (*)(mpfi_t,mpfi_t,mpfi_t)

void apply_binop(void (*binop)(mpfi_t,mpfi_t,mpfi_t),value x, value y, value z){
  location *loc_x, *loc_y, *loc_z;
  loc_x = value_to_location(x);
  loc_y = value_to_location(y);
  loc_z = value_to_location(z);
  mpfi_t *m_x, *m_y, *m_z;
  m_x = location_to_mpfi(loc_x);
  m_y = location_to_mpfi(loc_y);
  m_z = location_to_mpfi(loc_z);
  (*binop)(*m_z,*m_x,*m_y);
  return;
}

value add_(value x,value y){
  CAMLparam2(x,y);
  CAMLlocal1(z);
  z = new_value();
  apply_binop((BINOP) mpfi_add,x,y,z);
  CAMLreturn(z);
}

value sub_(value x,value y){
  CAMLparam2(x,y);
  CAMLlocal1(z);
  z = new_value();
  apply_binop((BINOP) mpfi_sub,x,y,z);
  CAMLreturn(z);
}

value mul_(value x,value y){
  CAMLparam2(x,y);
  CAMLlocal1(z);
  z = new_value();
  apply_binop((BINOP) mpfi_mul,x,y,z);
  CAMLreturn(z);
}

value div_(value x,value y){
  CAMLparam2(x,y);
  CAMLlocal1(z);
  z = new_value();
  apply_binop((BINOP) mpfi_div,x,y,z);
  CAMLreturn(z);
}

value union_(value x,value y){
  CAMLparam2(x,y);
  CAMLlocal1(z);
  z = new_value();
  apply_binop((BINOP) mpfi_union,x,y,z);
  CAMLreturn(z);
}

// mpfi -> mpfi -> int

#define MMOP int (*)(mpfi_t,mpfi_t)

int apply_mmop(int (*mmop)(mpfi_t,mpfi_t),value x, value y){
  location *loc_x,*loc_y;
  loc_x = value_to_location(x);
  loc_y = value_to_location(y);
  mpfi_t *m_x,*m_y;
  m_x = location_to_mpfi(loc_x);
  m_y = location_to_mpfi(loc_y);
  int dest = (*mmop)(*m_x,*m_y);
  return dest;
}

value cmp_(value x,value y){
  CAMLparam2(x,y);
  int dest = apply_mmop((MMOP) mpfi_cmp,x,y);
  CAMLreturn(Val_int(dest));
}

value is_inside_(value x,value y){
  CAMLparam2(x,y);
  int dest = apply_mmop((MMOP) mpfi_is_inside,x,y);
  CAMLreturn(Val_int(dest));
}
  
/* -------------------------------------------------------------------------- */
/*  Unary Operators                                                           */
/* -------------------------------------------------------------------------- */

// mpfi -> mpfi 

#define UNOP void (*)(mpfi_t,mpfi_t)

void apply_unop(void (*unop)(mpfi_t,mpfi_t),value x,value y){
  location *loc_x,*loc_y;
  loc_x = value_to_location(x);
  loc_y = value_to_location(y);
  mpfi_t *m_x,*m_y;
  m_x = location_to_mpfi(loc_x);
  m_y = location_to_mpfi(loc_y);
  (*unop)(*m_y,*m_x);
  return;
}

value neg_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_neg,x,y);
  CAMLreturn(y);
}

value inv_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_inv,x,y);
  CAMLreturn(y);
}

value sqr_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_sqr,x,y);
  CAMLreturn(y);
}

value sqrt_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_sqrt,x,y);
  CAMLreturn(y);
}

value abs_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_abs,x,y);
  CAMLreturn(y);
}

value log_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_log,x,y);
  CAMLreturn(y);
}

value sin_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_sin,x,y);
  CAMLreturn(y);
}

value cos_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_cos,x,y);
  CAMLreturn(y);
}

value tan_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_tan,x,y);
  CAMLreturn(y);
}

value asin_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_asin,x,y);
  CAMLreturn(y);
}

value acos_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_acos,x,y);
  CAMLreturn(y);
}

value atan_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unop((UNOP) mpfi_atan,x,y);
  CAMLreturn(y);
}

// mpfi -> mpfi (from float)

#define UNOPF void (*)(mpfr_t,mpfi_t)

void apply_unopf(void (*unopf)(mpfr_t,mpfi_t),value x,value y){
  location *loc_x,*loc_y;
  loc_x = value_to_location(x);
  loc_y = value_to_location(y);
  mpfi_t *m_x,*m_y;
  m_x = location_to_mpfi(loc_x);
  m_y = location_to_mpfi(loc_y);
  mpfr_t mpfr;
  mpfr_init(mpfr);
  (*unopf)(mpfr,*m_x);
  mpfi_set_fr(*m_y,mpfr);
  mpfr_clear(mpfr);  
  return;
}

value icenter_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unopf((UNOPF) mpfi_mid,x,y);
  CAMLreturn(y);
}

value iwidth_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unopf((UNOPF) mpfi_diam_abs,x,y);
  CAMLreturn(y);
}

value ilower_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unopf((UNOPF) mpfi_get_left,x,y);
  CAMLreturn(y);
}

value iupper_(value x){
  CAMLparam1(x);
  CAMLlocal1(y);
  y = new_value();
  apply_unopf((UNOPF) mpfi_get_right,x,y);
  CAMLreturn(y);
}

// mpfi -> mpfi * mpfi

value bisect_(value x){
  CAMLparam1(x);
  CAMLlocal3(tuple,y,z);
  y = new_value();
  z = new_value();
  tuple = caml_alloc_tuple(2);
  Store_field(tuple,0,y);
  Store_field(tuple,1,z);
  location *loc_x,*loc_y,*loc_z;
  loc_x = value_to_location(x);
  loc_y = value_to_location(y);
  loc_z = value_to_location(z);
  mpfi_t *m_x,*m_y,*m_z;
  m_x = location_to_mpfi(loc_x);
  m_y = location_to_mpfi(loc_y);
  m_z = location_to_mpfi(loc_z);
  mpfi_bisect(*m_y,*m_z,*m_x);
  CAMLreturn(tuple);
}

// mpfi -> bool

#define BUNOP int (*)(mpfi_t)

value apply_bool_unop(int (*bunop)(mpfi_t),value x){
  location* loc_x;
  loc_x = value_to_location(x);
  mpfi_t* m_x;
  m_x = location_to_mpfi(loc_x);
  int n = (*bunop)(*m_x);
  return (Val_bool(n));
}

value has_zero_(value x){
  CAMLparam1(x);
  CAMLreturn(apply_bool_unop((BUNOP) mpfi_has_zero,x));
}

// mpfi_is_zero doesn't work correctly!!!
/* value is_zero_(value x){ */
/*   CAMLparam1(x); */
/*   CAMLreturn(apply_bool_unop((BUNOP) mpfi_is_zero,x)); */
/* } */

value nan_p_(value x){
  CAMLparam1(x);
  CAMLreturn(apply_bool_unop((BUNOP) mpfi_nan_p,x));
}

value inf_p_(value x){
  CAMLparam1(x);
  CAMLreturn(apply_bool_unop((BUNOP) mpfi_inf_p,x));
}

value bounded_p_(value x){
  CAMLparam1(x);
  CAMLreturn(apply_bool_unop((BUNOP) mpfi_bounded_p,x));
}

value is_empty_(value x){
  CAMLparam1(x);
  CAMLreturn(apply_bool_unop((BUNOP) mpfi_is_empty,x));
}

value is_strictly_pos_(value x){
  CAMLparam1(x);
  CAMLreturn(apply_bool_unop((BUNOP) mpfi_is_strictly_pos,x));
}

// mpfi -> float

#define FUNOP void (*)(mpfr_t,mpfi_t)

value apply_float_unop(void (*funop)(mpfr_t,mpfi_t),value x){
  location *loc_x; 
  loc_x = value_to_location(x);
  mpfi_t* m_x;
  m_x = location_to_mpfi(loc_x);
  mpfr_t w;
  mpfr_init(w);
  (*funop)(w,*m_x);
  double d = mpfr_get_d(w,GMP_RNDN);
  mpfr_clear(w);
  return (caml_copy_double(d));
}

value width_(value x){
  CAMLparam1(x);
  CAMLreturn(apply_float_unop((FUNOP) mpfi_diam_abs,x));
}

value center_(value x){
  CAMLparam1(x);
  CAMLreturn(apply_float_unop((FUNOP) mpfi_mid,x));
}

value lower_(value x){
  CAMLparam1(x);
  CAMLreturn(apply_float_unop((FUNOP) mpfi_get_left,x));
}

value upper_(value x){
  CAMLparam1(x);
  CAMLreturn(apply_float_unop((FUNOP) mpfi_get_right,x));
}

// mpfi -> float -> bool

value mem_(value x, value d){
  CAMLparam2(x,d);
  double f = Double_val(d);
  location *loc_x;
  loc_x = value_to_location(x);
  mpfi_t* m;
  m = location_to_mpfi(loc_x);
  int n = mpfi_is_inside_d(f,*m);
  CAMLreturn(Val_bool(n));
}

/* -------------------------------------------------------------------------- */
/*  Constants                                                                 */
/* -------------------------------------------------------------------------- */

value pi_(value unit){
  CAMLparam1(unit);
  CAMLlocal1(pi);
  pi = new_value();
  location *loc_pi;
  loc_pi = value_to_location(pi);
  mpfi_t *m_pi;
  m_pi = location_to_mpfi(loc_pi);
  mpfi_const_pi(*m_pi);
  CAMLreturn(pi);
}
